library(stringr)
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.4     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ forcats 0.5.1
✓ readr   1.4.0     
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x readr::col_factor() masks scales::col_factor()
x purrr::discard()    masks scales::discard()
x tidyr::expand()     masks Matrix::expand()
x dplyr::filter()     masks stats::filter()
x dplyr::lag()        masks stats::lag()
x tidyr::pack()       masks Matrix::pack()
x tidyr::unpack()     masks Matrix::unpack()
response <- try(system('~/google-cloud-sdk/bin/gcloud projects list --quiet', intern = T))
projectid <- strsplit(response[2], " ")[[1]][1]
options(na.action = "na.fail") 
source("./dredge_functions.R")
dredge_and_subset <- function(data) {
  model <- lm(
    presence_ratio ~ foraging_niche + trophic_niche + is_noctural + pc1 + pc2 + pc3 + pc4,
    data=data
  )
  dredge_result <- dredge(model)
  summary(model.avg(dredge_result, subset = delta < 10))
}
load_species_data <- function(poolname) {
  filename <- str_replace('species_metrics_##POOL_NAME##.csv', '##POOL_NAME##', poolname)
  
  if (!file.exists(filename)) {
    column_name = str_replace_all('present_in_##POOL_NAME##_pool', '##POOL_NAME##', poolname)
    if (poolname == 'both') {
      column_name = 'present_in_both_pools'
    }
    
    sql <- str_replace_all("
         WITH species AS (
              SELECT 
                  scientific_name,
                  count(*) AS regional_pool_count,
                  countif(present_in_city = true) AS city_count
              FROM model.regional_species_filtered
              JOIN model.taxonomy USING (scientific_name)
              WHERE ##COL_SELECT## = true
              GROUP BY scientific_name
          )
          SELECT 
              species.*,
              body_morphspace.pc1,
              body_morphspace.pc2,
              body_morphspace.pc3,
              body_morphspace.pc4,
              trophic_niche,
              foraging_niche,
              is_noctural,
              ttc.cluster,
              ttc.cluster_half,
              ttc.cluster_double
          FROM species
          JOIN model.taxonomy USING (scientific_name)
          JOIN datasci.taxonomic_trait_clusters ttc USING (scientific_name)", '##COL_SELECT##', column_name)
  
    tb <- bq_project_query(projectid, sql)

    data <- bq_table_download(tb)
    write_csv(data, filename)
  }
  
  data <- read_csv(filename)
  
  data[is.na(data$foraging_niche),]$foraging_niche <- 'Omnivore'
  
  data$foraging_niche = as.factor(data$foraging_niche)
  data$trophic_niche = as.factor(data$trophic_niche)
  data$is_noctural = as.factor(data$is_noctural)
  data$presence_ratio = data$city_count / data$regional_pool_count
  
  data
}
merlin_data <- load_species_data('merlin')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)
merlin_data
merlin_result <- dredge_and_subset(merlin_data)
Fixed term is "(Intercept)"
merlin_summ <- model_summary('merlin', 'species', merlin_result)
merlin_summ
birdlife_data <- load_species_data('birdlife')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)
birdlife_result <- dredge_and_subset(birdlife_data)
Fixed term is "(Intercept)"
birdlife_summ <- model_summary('birdlife', 'species', birdlife_result)
birdlife_summ
both_data <- load_species_data('both')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)
both_result <- dredge_and_subset(both_data)
Fixed term is "(Intercept)"
both_summ <- model_summary('both', 'species', both_result)
both_summ
either_data <- load_species_data('either')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)
either_result <- dredge_and_subset(either_data)
Fixed term is "(Intercept)"
either_summ <- model_summary('either', 'species', either_result)
either_summ
Full result
all_species_results <- full_join(full_join(merlin_summ, birdlife_summ), full_join(both_summ, either_summ))
Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")
write_csv(all_species_results, "species_result.csv")

library(ggpubr)

jpeg(width = 1024, height = 1024)
myplot
dev.off()
null device 
          1 
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KbGlicmFyeShzdHJpbmdyKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKYGBgCgpgYGB7cn0KcmVzcG9uc2UgPC0gdHJ5KHN5c3RlbSgnfi9nb29nbGUtY2xvdWQtc2RrL2Jpbi9nY2xvdWQgcHJvamVjdHMgbGlzdCAtLXF1aWV0JywgaW50ZXJuID0gVCkpCnByb2plY3RpZCA8LSBzdHJzcGxpdChyZXNwb25zZVsyXSwgIiAiKVtbMV1dWzFdCmBgYAoKYGBge3J9Cm9wdGlvbnMobmEuYWN0aW9uID0gIm5hLmZhaWwiKSAKYGBgCgpgYGB7cn0Kc291cmNlKCIuL2RyZWRnZV9mdW5jdGlvbnMuUiIpCmBgYAoKYGBge3J9CmRyZWRnZV9hbmRfc3Vic2V0IDwtIGZ1bmN0aW9uKGRhdGEpIHsKICBtb2RlbCA8LSBsbSgKICAgIHByZXNlbmNlX3JhdGlvIH4gZm9yYWdpbmdfbmljaGUgKyB0cm9waGljX25pY2hlICsgaXNfbm9jdHVyYWwgKyBwYzEgKyBwYzIgKyBwYzMgKyBwYzQsCiAgICBkYXRhPWRhdGEKICApCiAgZHJlZGdlX3Jlc3VsdCA8LSBkcmVkZ2UobW9kZWwpCiAgc3VtbWFyeShtb2RlbC5hdmcoZHJlZGdlX3Jlc3VsdCwgc3Vic2V0ID0gZGVsdGEgPCAxMCkpCn0KYGBgCgpgYGB7cn0KbG9hZF9zcGVjaWVzX2RhdGEgPC0gZnVuY3Rpb24ocG9vbG5hbWUpIHsKICBmaWxlbmFtZSA8LSBzdHJfcmVwbGFjZSgnc3BlY2llc19tZXRyaWNzXyMjUE9PTF9OQU1FIyMuY3N2JywgJyMjUE9PTF9OQU1FIyMnLCBwb29sbmFtZSkKICAKICBpZiAoIWZpbGUuZXhpc3RzKGZpbGVuYW1lKSkgewogICAgY29sdW1uX25hbWUgPSBzdHJfcmVwbGFjZV9hbGwoJ3ByZXNlbnRfaW5fIyNQT09MX05BTUUjI19wb29sJywgJyMjUE9PTF9OQU1FIyMnLCBwb29sbmFtZSkKICAgIGlmIChwb29sbmFtZSA9PSAnYm90aCcpIHsKICAgICAgY29sdW1uX25hbWUgPSAncHJlc2VudF9pbl9ib3RoX3Bvb2xzJwogICAgfQogICAgCiAgICBzcWwgPC0gc3RyX3JlcGxhY2VfYWxsKCIKICAgICAgICAgV0lUSCBzcGVjaWVzIEFTICgKICAgICAgICAgICAgICBTRUxFQ1QgCiAgICAgICAgICAgICAgICAgIHNjaWVudGlmaWNfbmFtZSwKICAgICAgICAgICAgICAgICAgY291bnQoKikgQVMgcmVnaW9uYWxfcG9vbF9jb3VudCwKICAgICAgICAgICAgICAgICAgY291bnRpZihwcmVzZW50X2luX2NpdHkgPSB0cnVlKSBBUyBjaXR5X2NvdW50CiAgICAgICAgICAgICAgRlJPTSBtb2RlbC5yZWdpb25hbF9zcGVjaWVzX2ZpbHRlcmVkCiAgICAgICAgICAgICAgSk9JTiBtb2RlbC50YXhvbm9teSBVU0lORyAoc2NpZW50aWZpY19uYW1lKQogICAgICAgICAgICAgIFdIRVJFICMjQ09MX1NFTEVDVCMjID0gdHJ1ZQogICAgICAgICAgICAgIEdST1VQIEJZIHNjaWVudGlmaWNfbmFtZQogICAgICAgICAgKQogICAgICAgICAgU0VMRUNUIAogICAgICAgICAgICAgIHNwZWNpZXMuKiwKICAgICAgICAgICAgICBib2R5X21vcnBoc3BhY2UucGMxLAogICAgICAgICAgICAgIGJvZHlfbW9ycGhzcGFjZS5wYzIsCiAgICAgICAgICAgICAgYm9keV9tb3JwaHNwYWNlLnBjMywKICAgICAgICAgICAgICBib2R5X21vcnBoc3BhY2UucGM0LAogICAgICAgICAgICAgIHRyb3BoaWNfbmljaGUsCiAgICAgICAgICAgICAgZm9yYWdpbmdfbmljaGUsCiAgICAgICAgICAgICAgaXNfbm9jdHVyYWwsCiAgICAgICAgICAgICAgdHRjLmNsdXN0ZXIsCiAgICAgICAgICAgICAgdHRjLmNsdXN0ZXJfaGFsZiwKICAgICAgICAgICAgICB0dGMuY2x1c3Rlcl9kb3VibGUKICAgICAgICAgIEZST00gc3BlY2llcwogICAgICAgICAgSk9JTiBtb2RlbC50YXhvbm9teSBVU0lORyAoc2NpZW50aWZpY19uYW1lKQogICAgICAgICAgSk9JTiBkYXRhc2NpLnRheG9ub21pY190cmFpdF9jbHVzdGVycyB0dGMgVVNJTkcgKHNjaWVudGlmaWNfbmFtZSkiLCAnIyNDT0xfU0VMRUNUIyMnLCBjb2x1bW5fbmFtZSkKICAKICAgIHRiIDwtIGJxX3Byb2plY3RfcXVlcnkocHJvamVjdGlkLCBzcWwpCgogICAgZGF0YSA8LSBicV90YWJsZV9kb3dubG9hZCh0YikKICAgIHdyaXRlX2NzdihkYXRhLCBmaWxlbmFtZSkKICB9CiAgCiAgZGF0YSA8LSByZWFkX2NzdihmaWxlbmFtZSkKICAKICBkYXRhW2lzLm5hKGRhdGEkZm9yYWdpbmdfbmljaGUpLF0kZm9yYWdpbmdfbmljaGUgPC0gJ09tbml2b3JlJwogIAogIGRhdGEkZm9yYWdpbmdfbmljaGUgPSBhcy5mYWN0b3IoZGF0YSRmb3JhZ2luZ19uaWNoZSkKICBkYXRhJHRyb3BoaWNfbmljaGUgPSBhcy5mYWN0b3IoZGF0YSR0cm9waGljX25pY2hlKQogIGRhdGEkaXNfbm9jdHVyYWwgPSBhcy5mYWN0b3IoZGF0YSRpc19ub2N0dXJhbCkKICBkYXRhJHByZXNlbmNlX3JhdGlvID0gZGF0YSRjaXR5X2NvdW50IC8gZGF0YSRyZWdpb25hbF9wb29sX2NvdW50CiAgCiAgZGF0YQp9CmBgYAoKCmBgYHtyfQptZXJsaW5fZGF0YSA8LSBsb2FkX3NwZWNpZXNfZGF0YSgnbWVybGluJykKbWVybGluX2RhdGEKYGBgCgoKCmBgYHtyfQptZXJsaW5fcmVzdWx0IDwtIGRyZWRnZV9hbmRfc3Vic2V0KG1lcmxpbl9kYXRhKQptZXJsaW5fc3VtbSA8LSBtb2RlbF9zdW1tYXJ5KCdtZXJsaW4nLCAnc3BlY2llcycsIG1lcmxpbl9yZXN1bHQpCm1lcmxpbl9zdW1tCmBgYAoKYGBge3J9CmJpcmRsaWZlX2RhdGEgPC0gbG9hZF9zcGVjaWVzX2RhdGEoJ2JpcmRsaWZlJykKYmlyZGxpZmVfcmVzdWx0IDwtIGRyZWRnZV9hbmRfc3Vic2V0KGJpcmRsaWZlX2RhdGEpCmJpcmRsaWZlX3N1bW0gPC0gbW9kZWxfc3VtbWFyeSgnYmlyZGxpZmUnLCAnc3BlY2llcycsIGJpcmRsaWZlX3Jlc3VsdCkKYmlyZGxpZmVfc3VtbQpgYGAKCgpgYGB7cn0KYm90aF9kYXRhIDwtIGxvYWRfc3BlY2llc19kYXRhKCdib3RoJykKYm90aF9yZXN1bHQgPC0gZHJlZGdlX2FuZF9zdWJzZXQoYm90aF9kYXRhKQpib3RoX3N1bW0gPC0gbW9kZWxfc3VtbWFyeSgnYm90aCcsICdzcGVjaWVzJywgYm90aF9yZXN1bHQpCmJvdGhfc3VtbQpgYGAKCgpgYGB7cn0KZWl0aGVyX2RhdGEgPC0gbG9hZF9zcGVjaWVzX2RhdGEoJ2VpdGhlcicpCmVpdGhlcl9yZXN1bHQgPC0gZHJlZGdlX2FuZF9zdWJzZXQoZWl0aGVyX2RhdGEpCmVpdGhlcl9zdW1tIDwtIG1vZGVsX3N1bW1hcnkoJ2VpdGhlcicsICdzcGVjaWVzJywgZWl0aGVyX3Jlc3VsdCkKZWl0aGVyX3N1bW0KYGBgCgotLS0tLS0tLS0tLS0tLS0tLS0tLS0tCkZ1bGwgcmVzdWx0Ci0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KYGBge3J9CmFsbF9zcGVjaWVzX3Jlc3VsdHMgPC0gZnVsbF9qb2luKGZ1bGxfam9pbihtZXJsaW5fc3VtbSwgYmlyZGxpZmVfc3VtbSksIGZ1bGxfam9pbihib3RoX3N1bW0sIGVpdGhlcl9zdW1tKSkKd3JpdGVfY3N2KGFsbF9zcGVjaWVzX3Jlc3VsdHMsICJzcGVjaWVzX3Jlc3VsdC5jc3YiKQpgYGAKCmBgYHtyfQptZXJsaW5fZGF0YSRzb3VyY2UgPC0gJ01lcmxpbicKYmlyZGxpZmVfZGF0YSRzb3VyY2UgPC0gJ0JpcmRsaWZlJwpib3RoX2RhdGEkc291cmNlIDwtICdCb3RoJwplaXRoZXJfZGF0YSRzb3VyY2UgPC0gJ0VpdGhlcicKCnBsb3RfZGF0YSA8LSByYmluZChtZXJsaW5fZGF0YSwgYmlyZGxpZmVfZGF0YSwgYm90aF9kYXRhLCBlaXRoZXJfZGF0YSkKCmdncGxvdChwbG90X2RhdGEsIGFlcyh4ID0gcGMxLCB5ID0gcGMyLCBhbHBoYSA9IHByZXNlbmNlX3JhdGlvLCBjb2xvciA9IGFzLmZhY3Rvcihmb3JhZ2luZ19uaWNoZSkpKSArIGdlb21fcG9pbnQoKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKyBmYWNldF93cmFwKH4gc291cmNlKQpgYGAKYGBge3J9CmdncGxvdChwbG90X2RhdGEsIGFlcyh4ID0gcGMzLCB5ID0gcGM0LCBhbHBoYSA9IHByZXNlbmNlX3JhdGlvLCBjb2xvciA9IGFzLmZhY3Rvcihmb3JhZ2luZ19uaWNoZSkpKSArIGdlb21fcG9pbnQoKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKyBmYWNldF93cmFwKH4gc291cmNlKQpgYGAKYGBge3J9CmxpYnJhcnkoZ2dwdWJyKQpgYGAKCmBgYHtyfQptZXJsaW5fZGF0YSRuaWNoZSA8LSBwYXN0ZShtZXJsaW5fZGF0YSR0cm9waGljX25pY2hlLCAnICcsIG1lcmxpbl9kYXRhJGZvcmFnaW5nX25pY2hlLCAnICcsIG1lcmxpbl9kYXRhJGlzX25vY3R1cmFsLCAnICcsIG1lcmxpbl9kYXRhJGNsdXN0ZXIpCmBgYAoKYGBge3J9Cm15cGxvdCA8LSBnZ2FycmFuZ2UoCiAgbmNvbCA9IDIsCiAgbnJvdyA9IDIsCiAgZ2dwbG90KG1lcmxpbl9kYXRhLCBhZXMoeCA9IHBjMSwgeSA9IHBjMywgYWxwaGEgPSBwcmVzZW5jZV9yYXRpbywgY29sb3VyID0gbmljaGUpKSArIGdlb21fcG9pbnQoKSArIHRoZW1lX2J3KCkgKyAgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIiwgYXhpcy50aXRsZS54PWVsZW1lbnRfYmxhbmsoKSkgKyB5bGFiKCJTVCArIFBCIC0+IExUICsgU0IiKSwKICBnZ3Bsb3QobWVybGluX2RhdGEsIGFlcyh4ID0gcGMyLCB5ID0gcGMzLCBhbHBoYSA9IHByZXNlbmNlX3JhdGlvLCBjb2xvdXIgPSBuaWNoZSkpICsgZ2VvbV9wb2ludCgpICsgdGhlbWVfYncoKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbj0ibm9uZSIsIGF4aXMudGl0bGUueD1lbGVtZW50X2JsYW5rKCksIGF4aXMudGl0bGUueT1lbGVtZW50X2JsYW5rKCkpLAogIGdncGxvdChtZXJsaW5fZGF0YSwgYWVzKHggPSBwYzEsIHkgPSBwYzQsIGFscGhhID0gcHJlc2VuY2VfcmF0aW8sIGNvbG91ciA9IG5pY2hlKSkgKyBnZW9tX3BvaW50KCkgKyB0aGVtZV9idygpICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIikgKyB5bGFiKCJMVCArIFBCIC0+IFNUICsgU0IiKSArIHhsYWIoIkJvZHkgU2l6ZSIpLAogIGdncGxvdChtZXJsaW5fZGF0YSwgYWVzKHggPSBwYzIsIHkgPSBwYzQsIGFscGhhID0gcHJlc2VuY2VfcmF0aW8sIGNvbG91ciA9IG5pY2hlKSkgKyBnZW9tX3BvaW50KCkgKyB0aGVtZV9idygpICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIiwgYXhpcy50aXRsZS55PWVsZW1lbnRfYmxhbmsoKSkgKyB4bGFiKCJCZWFrIFNpemUiKQopCm15cGxvdApgYGAKYGBge3J9CmpwZWcod2lkdGggPSAxMDI0LCBoZWlnaHQgPSAxMDI0KQpteXBsb3QKZGV2Lm9mZigpCmBgYA==